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Abstract. Given a map / : Z — > Z and an initial argument a, we can iterate the 
map to get a finite set of iterates modulo a prime p. In particular, for a quadratic map 
f(z) = z 2 + c, c constant, work by Pollard suggests that this set should have length 
on the order of y/p. We give a heuristic argument that suggests that the statistical 
properties of this set might be very similar to the Birthday Problem random variable X n , 
for an n — p day year, and offer considerable experimental evidence that the limiting 
distribution of these set lengths, divided by y/p, for p < x as x — > oo, converges to the 
limiting distribution of X n /y/n, as n — s> oo. 



1. INTRODUCTION 

Let / € Z[z] be a polynomial and let a G Z. We define the orbit of a under / to be 

O f (a) = {f n (a) :n = 0,l,2,3,...}, 
and for each prime p we define the orbit modulo p of a under / to be 

O p f (a) = {f n (a) mod p : n = 0, 1, 2, 3, . . . }, 
where f n is the n th iterate of /, 

r=f°f° — °f, 

K v ' 

n 

and f°(a) = a. For a fixed / and a and a given prime p, let m p be the order of Oj(a). 

If / is a random map, i.e., a map chosen from the uniformly distributed set consisting of 
all maps from F p into F p (see Harris [6]), then the values of f n (a) are uniformly distributed 
for all n, and all a, and so the probability that /°(a),/ 1 (a),/ 2 (a), ...,/ fc (a) are all different 
is 

^ p — 1 p—2 p — k (p — 1)! 

p p p p k (p — k — 1)! ' 

since, once a is fixed, there are p—\ choices for J 1 (a), p — 2 choices for / 2 (a), and so on. 
Therefore, in this case the probability that (at least) two of /°(q!),/ 1 (q),/ 2 (q), . . . ,f k (a) 
are equal is 

(P)_! (P- 1 ) 1 
k p k (p — k — 1)! ' 

By an analogous argument, qi is also the probability that among k people, two people 
have the same birthday, where p is the number of days in a year. Framing this a little 
differently, we let the random variable X n be the number of times that we must sample 

l 
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(uniformly, with replacement) from the set {1,2,3,... , n} to get a repetition. Since it is 
known that the expected value of this variable is on the order of y/n, we look instead at 
the variable ^s. 

In light of the above heuristic, we might expect that for a fixed polynomial / and initial 
value a, ^£ will, on average, "behave" similarly to ^a. In particular, we might guess that 

the limiting distribution of for p < x, x — > oo, will be similar to the limiting distribution 
of — as n — > oo. We note that the above heuristic is not new; similar arguments have 
been given by Pollard [7], Bach [I], and Brent [3] to name a few, leading to conjectures 
that m p is on average ~ y^fp- 

We also consider a related question. For a fixed / £ Z[z], a E Z, let 

Q/,a( :r ) = {p < £ : f n {ot) = (mod p) for some n = 0, 1, 2, . . . }. 

That is, Qf t a(x) is the set of primes p less than or equal to x such that appears in the 
orbit modulo p of a under /. In particular, we are interested in the size of Qr a (x). Since, 
for a given prime p, the proportion of elements mod p in the orbit of a under / is 
we hypothesize that |Q/, a (x)| will grow at a rate proportional to Therefore, if we are 
correct that m p will grow at a rate proportional to y/p, we might expect that 



\QfA*)\ 



Em p yfx 
= C • \ : 
V log x 

p<x y to 



for some constant c € M. 

In the following we take an experimental approach to studying properties of the set 

For selected maps / and initial values a, we compute the orbits modulo p for all p < 2 25 . 
In particular, given these orbits we can find the moments of and the length of Qf )a (x). 
As we will demonstrate in the sections to follow, our results give strong support to the 
above heuristic, and lead us to make the following conjectures: 

Conjecture 1. Let f(z) = z 2 + c and a G Z be such that 

1) cGZ\{0,-2} 

2) a 7^±|(1 ±vT=ic), a ^±|(1 ±V-3-4c), a / 0, ±1 when c = — 1, 

and let the orbit length m p be as defined above. Then as x — >■ oo, the distribution of 

converges, independent of f and a, to a continuous distribution F(t) = 1 — e - * 2//2 ,i > 0. 
In particular, the r th moments of ^= are given by \i T = r(i — 2)(r — 4) • • • 2 for r even, and 

fj, r = r(r — 2)(r — 4) • • • 1 • yf% for r odd. 



The motivation for the result conjectured above is elaborated upon in section 2, and the 
need to include conditions (1) and (2) for both conjectures is explained in section 4. 
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Conjecture 2. Let f{z) = z 2 + c and a £ Z 6e suc/i £/ia£ conditions (1) and (2) o/ Conjec- 
ture 1 hold, and a 2 ^ —c. Define Qf >a (x) = {p < x : f n (a) = (mod p) for some n > 0}. 
TTien 

lim (\Q f>a (x)\^)=y/to 
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2. LENGTH OF THE ORBIT MODULO P AND THE BIRTHDAY PROBLEM 

Let Ek be the k th number drawn uniformly from the set {1, 2,3,... , n}, with replacement, 
and let X n be as defined in section 1. Then for k < n we have 

P{X n > k) = P{E\, . . . ,Ek all take different values) 



[ (1 - P(Ej = Ei for some i < j)) 

k / . -. \ rfe-1 

3 ~ 1 



i=2 



n 



exp 



J>g(l-j/ 



ii 



So as n — > oo, we have the following for < t < y/n: 

lim P{X n jxfn > t) = lim P(X n > ty/n) 

n— >oo n— >oo 

= lim exp 



lim exp 



log (1 - j/n) 

l<j<ty/ri 

l<j<ty/n k=l 



k 



where we have used the power series representation for log(l — j/n) in the third line. 
Switching the order of summation, and pulling the first term of the sum over k out of the 
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exponential, we have 



lim exp 



lim exp 



E H n 

l<j<ty/n 



2n 



lim exp 

n— >oo 



lim exp 



E E 



U/n ) 
k 



exp 



V lim O 

— J n— >oo 



k=l 



fk+2 



kn k / 2 



k=2 l<j<t,J7i 

~Y [^kj2 

-t 2 /2 



where the second line follows because, in general, J2i<j<mJ k 1S a polynomial in m of degree 
k + 1, and the third line, where we have brought the limit inside the sum, follows from 
the Monotone Convergence Theorem. Therefore lim n _ 5 . 00 P{X n /^/n < t) = 1 — e~ l / 2 , so 
we see that the distribution of ^= converges to a distribution function F(t) = 1 — e _< / 2 , 

which has an associated density function f(t) = F'(t) = te~ t2 1 2 . To support our conjecture 
in section 1 — that F(t) is the limiting distribution of as x — > 00 — we compare the 

moments of which we compute in section 5 for large x, to the limiting moments of 

as n — > 00. With the limiting density function f(t) of ^= in hand we can derive a general 

expression for the r th moment: 



/•oo 

/ t r f(t) 

JO 



dt = I f +1 e~ t2 / 2 dt 




r -t 2 /2 



-t r e 



oo roc 

+ r / f "V* 2 / 2 dt 
Jo 



■ / f-'e- 12 / 2 dt, 
'0 



where r applications of L'Hospital's rule give us for the — t r e * 2//2 term. We continue 
applying integration by parts as above until we get 



r is even, 



poo 

fx r = r(r-2)(r-4)---2- / te~ t2/2 dt if 

Jo 

poo 

(r-2)(r-4)---l- / e~* 2/2 dt if r is odd, 
Jo 
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The first integral above evaluates to — e~ l > 2 = 1, and the second integral we evaluate 
as follows: 

00 

2 , 



I = / e- t2 l 2 dt 
Jo 

(2/) 2 = e -* 2 /2 dt V 



00 /■oo 



e~ x /z dx- / e~ y /z dy 

■00 J —00 

00 poo 

/ e-^+^dxdy 
■00 J —00 

OO /"27T 

/ re" r2/2 drd# = 2vr 

r=0 ^9=0 
7T 

Therefore the r th moments of the limiting distribution of as n — >• 00, are given by 

\x r = r(r — 2)(r — 4) • • • 2 if r is even, 

H r = r(r - 2){r - A) ■ ■ ■ I ■ if r is odd. 



For the first four moments this gives us fii = yvr/2, fi2 = 2, {13 = 2>^/tt/2, \i± = 8. 
Therefore, to support our claim in Conjecture 1 we must provide evidence that the moments 
of ^= are converging, as x — > 00, to the moments fj, r above. In our computations we use 

the following expression for the r th moments of 

|{p<x}|^VVp 



3. ITERATES OF / CONGRUENT TO ZERO MODULO P 

In this section we consider the quantity |S/ iQ (^)|-^^, as defined in Section 1. Assuming 
that the probability that G Oj(a) is and that M\ will converge to y/ir/2, we define 

G(x) = I] P <x ~^T> an ^ mak e a guess that 



(1) lim (iQfA*)^) 



lim G(x). 
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If we let tt(x) = X^fc<x a (^)' where a(k) = 1 if k is prime and otherwise, and define 
/(x) = then Stieltjes- integration by parts gives 

i3/2 ' 

logx^ 1 Um /logs / ^g) d 



(2) =► hm Hep^J- 

x— >oc 1 ./nr. — / . fpi 



x^oo I ,/x x fp j x-too 1 2Jx I i 3 / 2 



Since 7r(x) is bounded by the inequality 

x 



< 7T X < 



log x + 2 log x — 4 ' 

for x > 55 [9], we have 

lim I * [ (-^ )dt) < lim / ^§dt 

*-**> \ V| jf Vv / i(logt + 2)/ I ^ool2^ jf t 3 / 2 

< lim I — I -= }dt 



M / Vv^(iogt-4)y 

\ log x J 2 

The bounding limits above have the indeterminate form — , since n „ 1 -■ > n „ 1 — — > 

° oo' Vt(logt-4) V*(logt + 2) 

j for i > 55 and 1/tdt diverges as x — > oo. Therefore we can apply L'Hopital's Rule to 
the bounding limits, giving us 

l=hm P|^-Ulim f_i r-#^"l<lim fe/ dt 

x-s>oo \log z x - 4/ y-y/x^ogx + 2) log —2 / z-+°o I 2vx J t 6 ' 2 

< Hm f-=^ ^log 2 x\ = Hm / log 2 x \ 

x^oo y^/x^logx - 4) log -2 y a;-+tx5 ylog^x - 6 log a; + 8J 

Therefore lim^oo G(x) = 2 • \Jtt/2 = \/2n, so our guess ([1]) becomes 

(3) lim (\Q f ( x )\^) =V2^. 

As we test our hypothesis, it should be kept in mind that lima^oo G{x) converges very 
slowly. Since the x values for which | Qf 7 a (x) \ -^S- can actually be computed (in a reasonable 
amount of time) are relatively small, the largest being 2 27 , we compare our computations 
to G(x), rather than the limit y/2~7r. 
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4. SOME SPECIAL CASES 

In this paper we consider polynomials of the form f(z) = z 2 + c, with z, c G Z, and initial 
argument values a G Z. But for certain /, a pairs we find that we end up with a finite (over 
Z) orbit, a condition which is clearly incompatible with our hypotheses outlined in sections 
2 and 3, since m p will have a fixed bound for all primes p. In this section we classify these 
exceptional pairs /, a. 

Proposition 1. Let Of (a) = {f n (a) : n = 0, 1, 2, 3, . . . } be the orbit of a under f , where 
f(z) = z 2 + c, c G Z, and a G Z. T/ien Of (a) is finite if and only if one of the following 
hold: 

i) a =±|(1 ±v / I = 45) 
ii) a =±|(1 ±V-3-4c) 
m) a G {0,1, -1} and c e {0, -1,-2}. 

Proof. First we prove the converse, which is easier. Assumption (i) gives us the solutions 
to a 2 ± a + c = 0, and this equation implies that a 2 + c = dba, which implies that the 
orbit is finite. Assumption (ii) gives the solutions to a 2 ± a + c + 1 = 0, and this equation 
implies that a 2 + c = ±a — 1. With one more iteration we get 

(a 2 + c) 2 + c = (±a - l) 2 + c = a 2 =F 2a + 1 - a 2 ± a - 1 = =f2q ± a = ±a, 

which again implies that the orbit is finite. As for (Hi), testing all possible a, c combina- 
tions will quickly convince the reader that the orbits are finite in all cases. 

Now suppose that Of (a) is finite. First we make some simplifications. Since the orbits 
of a and —a will be identical except for the sign of the first element / = a, we may 
consider only non-negative values of a. Also, since it is obvious that c G {0, —1} will have 
infinite orbit for a > 2, and that c > 1 will have infinite orbit for all a, we consider only 
c < —2. We claim that Of (a) finite implies \[—c — 1 < a < \/—c+ 1. If this were not true, 
then we would have either a = [\/— c] + b or a = [\^—c\ — b for some b G N, giving us 

a = [V=cl +b a 2 + c = ( [v 7 ^! ) 2 + 2fr[ v /Z cl + b 2 + c> [x/^c] + b 

a=[yf^\-b a 2 + c= (Lv /Z cJ) 2 -2&Lv /Z cJ +^ 2 + c< -26[v /Z cJ +c. 

The first of these immediately implies that the iterates of / are unbounded since they are 
strictly increasing. In the second case iterating once more gives us 

(a 2 + c) 2 + c> 46 2 Lv^J 2 - 46cLv /Z cJ + c 2 > Lv^J + b 

where the inequality reverses since a 2 + c < — 26 [\/— cj + c < 0, and the second inequality 
follows since c < —2. Again we can conclude that the iterates of / are unbounded, and 
so we have shown that Of (a) finite implies y/^c — 1 < a < \f—c + 1. For any c, there 
are at most two integers that satisfy the preceding inequality, Lv^^cJ an d \\/~ c| , so any 
member of Of (a) must be one of cj , ± [V — c l > since otherwise the iterates of / will 

be unbounded. Since we know a G Of (a), the condition above implies that Of (a) C 
{a, —a, a — 1, —a — 1} or Of (a) C {a, —a, a + 1, —a + 1} . However, we can rule out the 
latter case since 
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a + c = ±a> + 1 
(a 2 + c) 2 + c = ±3a + 2 
((a 2 + c) 2 + c) 2 + c = 7a 2 ± 13a + 5 > 2a + 5 > ±a + 1 > ±a, 

where the first inequality follows since in this case c < — 2 ==>■ a > 2. Therefore the 
iterates are unbounded in this case, and we are left with the following: 

a 2 + c = ±a — 1 
a 2 ±a + c+ l = 
a =±|(1 ±V-3-4c). 

□ 

This proposition is the basis for the second condition necessary for Conjectures 1 and 2; 
we now turn to the first condition, that c ^ {0, —2}. For c = it is immediately clear that 
|Q/,a(^)| will not grow as expected, since we'll have p £ Qf, a {x) if and only if jd divides 
a. On the other hand, the length m p of the orbit modulo p will grow much faster than 
we expect, for both c = and c = —2. Vasiga and Shallit [TT] studied these cases, showing 
that, for a given prime p, if (p — l)/2 is prime and 2 is a primitive root modulo (p — l)/2, 
then ^o<«<p m p * s a * l eas t on the order of p 2 . Heuristics by Hardy and Littlewood [5], 
along with Artin's conjecture, suggest that the number of primes less than x that satisfy 
this property is on the order of x/(logx) 2 , and thus the density of these primes is on the 
order of 1/ log x. If we sum p 2 , for p < x, and multiply by 1/ log x, we get something on the 
order of x 3 /(logx) 2 , and dividing this by the sum ^2 p<x J2o<a<p ^ ~ x 2 /logx gives us an 
average orbit length of ~ x/logx. Note that this estimate only takes into account primes 
with the aforementioned property, and assumes that all other primes have orbit length 0, 
so we should expect this to be a low estimate. Indeed, the limited experimentation we did 
on this question suggests that the average orbit length is closer to x/(logx) 3 / 4 . 

Finally, Conjecture 2 requires an additional condition, that a 2 ^ —c. If we disregard this 
condition then we will have cases where G Qf, a (x) for all p, which is clearly incongruent 
with our claim. To see that the / = a is the only iterate whose square can be equal to 
— c, suppose that the contrary were true, i.e. that we have (f 1 ) 2 = —c for some I € Z, then, 
letting f k = we have 

((/ fc ) 2 + c) 2 + c = 
(/ fc ) 4 + 2c(/ fc ) 2 + c 2 + c = 
c 2 + (2(/ fe ) 2 + l)c+(/ fc ) 4 = 

Therefore the quadratic formula gives us 




a + c = ±a or 



a±a + c = or 



a =±i(l ±V^^c) or 



c = 



-2(/ fc ) 2 - 1 ± s/Wf + 1 



2 
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which is not an integer unless f k = 0, in which case (f 1 ) 2 = c 2 = — c ==>■ c = — 1. It is easy 
to see that this implies a E {0, 1}, and this case has already been excluded by Proposition 
l(iii). 



5. RESULTS 

First we consider the first, second, third and fourth moments of as discussed in 

VP 

section 2, for f(z) = z 2 + c, where c = ±1, +2, ±3, and initial arguments a = 1, 2, . . . , 9. 
Of these we can exclude a = 1,2 when f(z) = z 2 — 3, and a = 1 when f(z) = z 2 — 1, 
because these (/,a) combinations have finite orbits, as discussed above. For the other 42 
combinations, we find that our experimental results support our hypotheses very well. For 
the first moment we expected the limit to be y ir/2 = 1.25331413..., and for all (f,a) 
tested, M\ was between 1.25138 and 1.25351 for x = 2 25 , with an average value of 1.25279. 
Table [1] gives these figures along with the standard deviation of the set of results for each 
moment. It also shows the mean, standard deviation, minimum, and maximum of the 
set {\^j2-M-i\ : x = 2 25 , for (/, a) tested}, and similarly for the second, third and fourth 
moments. Our complete results are depicted graphically in Figures 1-4, for the first, second, 
third, and fourth moments, respectively. In each of these graphs the heavier red curve is 
the respective moment of X n /^/n, for n = x. Notice that the y-axes of these graphs are 
not scaled equally with respect to each other (they are stretched by a factor of two for 
each subsequent moment graph), so if we're interested in comparing how quickly two of 
the moments converge, Table [T] will be more helpful. 



Table 1. Moments of "^f= for x = 2 25 and distance from 
predicted limit. 





mean 


stand dev 


min 


max 




1.252795789 


0.000518158 


1.251387582 


1.253505370 


y / V2-A/i 


0.000544827 


0.000490241 


0.000000052 


0.001926555 


M 2 


1.998325027 


0.001690776 


1.993860194 


2.000539507 


2-M 2 


0.001810403 


0.001544894 


0.000034079 


0.006139806 


M 3 


3.755044605 


0.004998323 


3.742341997 


3.762285912 


\^^JV/2-M■ 3 \ 


0.005419269 


0.004427558 


0.000121838 


0.017600415 


M 4 


7.985456401 


0.014915109 


7.948531018 


8.008817811 


8-M 4 


0.016278430 


0.012999594 


0.000149649 


0.051468982 



75 ~ 1.25331413731550 
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Figure 5. Histogram, H(t), of the distribution of ^| (blue) for x = 10 8 , 
f(z) = z 2 + 1, a = 3, superimposed on the graph of f(t) = te~^ I 2 (red). 
H(t :wk<t< w{k + 1)) = l{p - 1 ° 8: "g ( { ^ | ) | < " (fc+1)}l , k e N. Each bar 
of the histogram has width w ~ 5.6/800. 




t 

The apparent common limit of the moments of ^= and ^= suggests that the limiting 
distributions of as x — > oo, and the random variable X n /y/n, as n — > oo, are the 
same. For the variable -4= we showed in section 2 that as n — > oo, the distribution 

P(X n /y/n < t) converges to the function F(t) = 1 — e~* I 2 . As the histogram in Figure 
[5] shows, the density function f(t) = F'(t) approximates quite well the distribution of 

for x = 10 8 , f(z) = z 2 + 1, a = 3. These results give considerable support to our first 
conjecture, stated in section 1. 

To test the hypothesis discussed in section 3, we compute | Q/, Q (x)| ^j=r, for (f,a) as 

described above and x S {2, 2 2 , . . . , 2 27 }. Table [2] shows that, although our results are still 
fairly widely dispersed at x = 2 27 , the average of the results for this x value is very close to 
G(x), and the standard deviation is decreasing in general as x increases, as is the error of 
the mean from G(x). As we mentioned earlier, huL^oo G(x) converges very slowly, and as 
the table shows, even for x as large as 2 27 we still have \G(x) — V2tt\ ~ 0.36, so we are not 
too surprised to see such a wide range in our results for this x value. That is, intuitively, 
it seems we should not expect our results to be very tightly grouped until we are close to 
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the limiting value, \^2tt. Figure [6] gives a graphical representation of all (/, a) tested, for 
x from 4 to 2 27 . On this graph the red and blue lines are G(x) and the mean from Table 
[2j respectively. From this data, it seems reasonable to suppose that \Qf, a {x)\^§ will 

eventually converge to \^2tt, independent of /, a, and so we make our second conjecture 
as stated in section 1. 

In his paper Variation of Periods Modulo p in Arithmetic Dynamics, Silverman carries 
out computations that lead to a conjecture (in a more general setting) that under certain 
restrictions the set {p : m p < p 1 / 2-6 } will have density for e > [10]. This conjecture 
agrees with our own results, and in fact if Conjecture 1 were proven, a less general version 
of Silverman's conjecture would readily follow. Computations of a similar nature to ours 
were also carried out in Periods of Rational Maps Modulo Primes by Benedetto, et al, with 
results that are compatible with our own [2]. 



Figure 6. Graphs of Qj a (x) l ^=- for all 42 (f,a) combinations tested 

(thinner lines), the mean of these graphs (thick blue line), and our guess 
G(x) (thick red line). 




9 10 11 12 13 14 15 16 17 16 

log 2 a; 



19 20 21 22 23 24 25 26 27 
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Table 2. A comparison of our experimental results to our guess, G(x). 



|S/,a(*)|^ for all (J, a) tested 



X 


Lr[X ) 


mean 


stand dev 


minimum 


maximum 


— mean) 


2 io 


3.46925 


3.00157 


0.51687 


1.73287 


4.11556 


0.46767 


2 11 


3.45003 


2.87221 


0.46933 


2.02178 


3.70660 


0.57781 


2 12 


3.42304 


2.79734 


0.34646 


2.07944 


3.37909 


0.62570 


nl3 


3.37313 


2.83502 


0.37519 


1.79203 


3.68363 


0.53811 


2 14 


3.32854 


2.87187 


0.30915 


1.89532 


3.56321 


0.45667 


2 15 


3.27415 


2.93202 


0.35071 


2.01030 


3.44622 


0.34213 


2 16 


3.22425 


2.97785 


0.33397 


2.20941 


3.63902 


0.24640 


2 


3.17737 


3.03080 


0.28861 


2.44107 


3.48260 


0.14657 


2 18 


3.13140 


3.04548 


0.18849 


2.55869 


3.38722 


0.08593 


2 19 


3.09045 


3.01797 


0.21004 


2.54637 


3.32847 


0.07248 


2 20 


3.05137 


2.93775 


0.17602 


2.63992 


3.27620 


0.11362 


2 21 


3.01654 


2.92737 


0.15786 


2.65359 


3.28683 


0.08917 


2 22 


2.98475 


2.94397 


0.12988 


2.71031 


3.21664 


0.04077 


2 23 


2.95589 


2.91037 


0.09402 


2.72467 


3.11548 


0.04552 


2 24 


2.92986 


2.88060 
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2.76176 


3.07449 


0.04925 


2 25 


2.90646 


2.88289 


0.05445 


2.77612 


3.02741 
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2 26 


2.88509 


2.87751 


0.05869 


2.71911 


3.01390 


0.00759 


2 27 


2.86578 


2.86821 


0.05418 


2.74621 


3.00790 


0.00243 


a log a- 


y/ir/2 
<x ^/p 
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